subroutine gl_moves
  use rbf_type
  use variable , only : nn,nodes,q
  implicit none
  real(dp)    , dimension(:,:) ,allocatable  :: q_new
  allocate (q_new(5,nodes))

  

  call bulk_moves (q_new)
  call surf_moves (q_new)
!  call surface_nano(a_new)

!  q(:,nn+1:nodes) = q_new(:,nn+1:nodes)
  q(:,1:nodes) = q_new(:,1:nodes)

  deallocate (q_new)
end subroutine gl_moves

subroutine bulk_moves(q_new)

  use rbf_type
  use variable
  use random_number2
  use ranftest
  use ifport
  use patch_module, only: inside_nano,surface_nano
  implicit none

  integer(i4b)           :: i,imc,nclus,k0,k1,k2,ii,k,j,l
  real(dp)    , dimension(5)  :: nrandom,a_bk,q_tmp

  real(dp)               :: du,invtemp,scalar,coef,gamma
  real(dp)    , dimension(3,3,6) :: q_full
  
  real(dp)    , dimension(3,3)   ::  q_imc,q_tmp_new

!----------------------------------------------------------------------------------  
  real(dp)    , dimension(5,nodes)   :: q_new
  real(dp)  ::  Elx,Ely,Elz,Eltot,LdG,delta,doubledot,dotq,dotd,Elchi,levi,levi1,levi2,dq,grad,A5,deltad,anch_deg,partial
  real(dp)  ::  d2Q,doubled,partialcross,El_l2_l4



  do imc =nn+1,nodes


     nanosurf:if(.not.surface_nano(imc).and..not.(inside_nano(imc)))then
        
        call q_qtensor (q(:,imc),q_imc(:,:))
        
        do i=1,6
           call q_qtensor (q(:,ix1(i,imc)),q_full (:,:,i)) 
        enddo
        
        doubled=doubledot(q_imc(:,:))
     
        ii = 0
     
        do i=1,3
           do j=i,3
          
              ii = ii + 1
              if(i.eq.j) then
                 delta = 1.d0/3.d0
              else
                 delta =0.d0
              endif
        
           
              dotd = dotq (q_imc(:,:),ii) 

              LdG = (1.d0 - UldG/3.d0)*q_imc (i,j) - ULdG * dotd  &
                   + Uldg* doubled * (q_imc(i,j) + delta )

              Elx = ddx2 *(q_full (i,j,1) - 2.d0*q_imc(i,j) + q_full (i,j,2) )
              Ely = ddy2 *(q_full (i,j,3) - 2.d0*q_imc(i,j) + q_full (i,j,4) )
              Elz = ddz2 *(q_full (i,j,5) - 2.d0*q_imc(i,j) + q_full (i,j,6) )

              Eltot = Elx + Ely + Elz

              if(lhat(2).ne.0.d0)then
                 if (ii.eq.1)then
                    d2Q= 2.d0*(partial(1,1,1,1,imc,q_imc(1,1)) + partial (1,2,1,2,imc,0.d0,ii)+partial (1,3,1,3,imc,0.d0))

                 else if (ii.eq.2)then
                    d2Q=   partial (2,1,1,1,imc,0.d0)       + partial (2,2,1,2,imc,q_imc(1,2))  &
                         + partial (2,3,1,3,imc,0.d0)       + partial (1,1,2,1,imc,q_imc(2,1))  &
                         + partial (1,2,2,2,imc,0.d0)       + partial (1,3,2,3,imc,0.d0)
                 else if ( ii.eq.3)then
                    d2Q=   partial (3,1,1,1,imc,0.d0)       + partial (3,2,1,2,imc,0.d0)        &
                         + partial (3,3,1,3,imc,q_imc(1,3)) + partial (1,1,3,1,imc,q_imc(3,1))  &
                         + partial (1,2,3,2,imc,0.d0)       + partial (1,3,3,3,imc,0.d0)
                 else if ( ii.eq.4)then
                    d2Q= 2.d0*(partial(2,1,2,1,imc,0.d0)    + partial (2,2,2,2,imc,q_imc(2,2))+partial (2,3,2,3,imc,0.d0))
                 else if ( ii.eq.5)then
                    d2Q=   partial (3,1,2,1,imc,0.d0)       + partial (3,2,2,2,imc,0.d0)  &
                         + partial (3,3,2,3,imc,q_imc(2,3)) + partial (2,1,3,1,imc,0.d0)  &
                         + partial (2,2,3,2,imc,q_imc(3,2)) + partial (2,3,3,3,imc,0.d0)
                 else if ( ii.eq.6)then
                    d2Q= 2.d0*(partial(3,1,3,1,imc,0.d0)    + partial (3,2,3,2,imc,0.d0)   +partial (3,3,3,3,imc,q_imc(3,3)))
                 endif

                 if (delta.eq.0.d0)then
                    partialcross =0.d0
                 else
                    partialcross =delta * ( &
                         partial(1,1,1,1,imc,q_imc(1,1)) + partial (2,2,2,2,imc,q_imc(2,2))+partial (3,3,3,3,imc,q_imc(3,3)) + &
                         partial(1,2,1,2,imc,0.d0) + partial (1,3,1,3,imc,0.d0)+partial (2,1,2,1,imc,0.d0) + &
                         partial(2,3,2,3,imc,0.d0) + partial (3,1,3,1,imc,0.d0)+partial (3,2,3,2,imc,0.d0) )
                 endif

                 El_l2_l4 = (lhat(2) + lhat(4))* (0.5d0*d2Q -  partialcross)
              else
                 El_l2_l4 = 0.d0
              endif

       !    full = cx * (Qold[i+1][j][k][0] - 2*Qold[i][j][k][0] + Qold[i-1][j][k][0]) + \
       !           cy * (Qold[i][j+1][k][3] - 2*Qold[i][j][k][3] + Qold[i][j-1][k][3]) + \
       !           cz * (Qold[i][j][k+1][5] - 2*Qold[i][j][k][5] + Qold[i][j][k-1][5]) + \
       !        2*(idxdy*(Qold[i+1][j+1][k][1] - Qold[i+1][j-1][k][1] - Qold[i-1][j+1][k][1] + Qold[i-1][j-1][k][1]) +\
       !           idxdz*(Qold[i+1][j][k+1][2] - Qold[i+1][j][k-1][2] - Qold[i-1][j][k+1][2] + Qold[i-1][j][k-1][2]) +\
       !           idydz*(Qold[i][j+1][k+1][4] - Qold[i][j+1][k-1][4] - Qold[i][j-1][k+1][4] + Qold[i][j-1][k-1][4]) );

        
       !       Elchi = 0.d0       
       !       El_l2_l4 = (lhat(2) + lhat(4))* (0.5d0*d2Q -  partialcross)
!              - dt*( LdG[n] - Elast[n] - (L2+L3)*(0.5*d2Q[n] - delta3(n)*full) );
              q_tmp_new (i,j)=  q_imc (i,j) + dmc(1) * ( -LdG  + Eltot + El_l2_l4  )
      
           enddo
        enddo
     !   write(*,*)q_tmp_new (1,1) + q_tmp_new (2,2) + q_tmp_new (3,3)
        
    call qtensor_q (q_tmp_new,q_new(:,imc))
    
   else
   q_new(:,imc)  = q(:,imc)

endif nanosurf
enddo

return
 456 format (I6,5f12.8)

end subroutine bulk_moves



function partial(k,l,i,j,imc,q_imc)
use rbf_type
use variable, only :ddx2,ddy2,ddz2,q,idxdy,idxdz,idydz
use variable, only :ix1,ix2,ix3,ix4,ix5
implicit none
integer i,j,k,l,imc,ii
real(dp) :: partial
real(dp)    , dimension(3,3,4) :: q_full
!real(dp)    :: q_1,q_2
real(dp)     :: q_1,q_2,q_3,q_4
 real(dp)      ::  q_imc
q_full = 0.d0
if (i.eq.1.and.j.eq.1) ii = 1 
if (i.eq.1.and.j.eq.2) ii = 2
if (i.eq.1.and.j.eq.3) ii = 3

if (i.eq.2.and.j.eq.2) ii = 4
if (i.eq.2.and.j.eq.3) ii = 5
if (i.eq.3.and.j.eq.3) ii = 6

if (i.eq.2.and.j.eq.1) ii = 2
if (i.eq.3.and.j.eq.1) ii = 3
if (i.eq.3.and.j.eq.2) ii = 5

if (k.eq.1.and.l.eq.1)then
   if(ii.ne.6)then
      partial  = ddx2 *(q(ii,ix1(1,imc)) - 2.d0*q_imc + q(ii,ix1(2,imc)) )
   else
      partial  = ddx2 *( (-q(1,ix1(1,imc))-q(4,ix1(1,imc)) ) - 2.d0*q_imc + (-q(1,ix1(2,imc))-q(4,ix1(2,imc)) ))
   endif
else if (k.eq.2.and.l.eq.2)then
   if(ii.ne.6)then
      partial  = ddy2 *(q(ii,ix1(3,imc)) - 2.d0*q_imc + q(ii,ix1(4,imc)) )
   else
      partial  = ddy2 *( (-q(1,ix1(3,imc))-q(4,ix1(3,imc)) ) - 2.d0*q_imc + (-q(1,ix1(4,imc))-q(4,ix1(4,imc)) ))
   endif

else if (k.eq.3.and.l.eq.3)then
  if(ii.ne.6)then
      partial  = ddz2 *(q(ii,ix1(5,imc)) - 2.d0*q_imc + q(ii,ix1(6,imc)) )
   else
      partial  = ddz2 *( (-q(1,ix1(5,imc))-q(4,ix1(5,imc)) ) - 2.d0*q_imc + (-q(1,ix1(6,imc))-q(4,ix1(6,imc)) ))
   endif
endif


if((k.eq.1.and.l.eq.2).or.(k.eq.2.and.l.eq.1))then

!   call q_qtensor (q(:,ix3(2,imc)),q_full (:,:,1))  ! i+1,j+1
!   call q_qtensor (q(:,ix2(2,imc)),q_full (:,:,2))  ! i+1,j-1
!   call q_qtensor (q(:,ix3(1,imc)),q_full (:,:,3))  ! i-1,j+1
!   call q_qtensor (q(:,ix2(1,imc)),q_full (:,:,4))  ! i-1,j-1

   if(ii.ne.6)then
      q_1 = q(ii,ix3(2,imc))     ! i+1,j+1
      q_2 = q(ii,ix2(2,imc))     ! i+1,j-1
      q_3 = q(ii,ix3(1,imc))     ! i-1,j+1
      q_4 = q(ii,ix2(1,imc))     ! i-1,j-1
   else
      q_1 = -q(1,ix3(2,imc)) -q(4,ix3(2,imc))    ! i+1,j+1
      q_2 = -q(1,ix2(2,imc)) -q(4,ix2(2,imc))    ! i+1,j-1
      q_3 = -q(1,ix3(1,imc)) -q(4,ix3(1,imc))    ! i-1,j+1
      q_4 = -q(1,ix2(1,imc)) -q(4,ix2(1,imc))    ! i-1,j-1
   endif

   partial = idxdy *(q_1 - q_2 - q_3  + q_4 )

 !  partial = idxdy *(q_full (i,j,1) - q_full (i,j,2) - q_full (i,j,3) + q_full (i,j,4))

!idxdy*(Qold[i+1][j+1][k][1] - Qold[i+1][j-1][k][1] - Qold[i-1][j+1][k][1] + Qold[i-1][j-1][k][1]) 
else if((k.eq.1.and.l.eq.3).or.(k.eq.3.and.l.eq.1))then

   if(ii.ne.6)then
      q_1 = q(ii,ix5(2,imc))     ! i+1,j+1
      q_2 = q(ii,ix4(2,imc))     ! i+1,j-1
      q_3 = q(ii,ix5(1,imc))     ! i-1,j+1
      q_4 = q(ii,ix4(1,imc))     ! i-1,j-1
   else
      q_1 = -q(1,ix5(2,imc)) -q(4,ix5(2,imc))    ! i+1,j+1
      q_2 = -q(1,ix4(2,imc)) -q(4,ix4(2,imc))    ! i+1,j-1
      q_3 = -q(1,ix5(1,imc)) -q(4,ix5(1,imc))    ! i-1,j+1
      q_4 = -q(1,ix4(1,imc)) -q(4,ix4(1,imc))    ! i-1,j-1
   endif

   partial = idxdz *(q_1 - q_2 - q_3  + q_4 )

!   call q_qtensor (q(:,ix5(2,imc)),q_full (:,:,1))  ! i+1,k+1
!   call q_qtensor (q(:,ix4(2,imc)),q_full (:,:,2))  ! i+1,k-1
!   call q_qtensor (q(:,ix5(1,imc)),q_full (:,:,3))  ! i-1,k+1
!   call q_qtensor (q(:,ix4(1,imc)),q_full (:,:,4))  ! i-1,k-1

!   partial = idxdz *(q_full (i,j,1) - q_full (i,j,2) - q_full (i,j,3) + q_full (i,j,4))
!idxdz*(Qold[i+1][j][k+1][2] - Qold[i+1][j][k-1][2] - Qold[i-1][j][k+1][2] + Qold[i-1][j][k-1][2]))
else if((k.eq.2.and.l.eq.3).or.(k.eq.3.and.l.eq.2))then

   if(ii.ne.6)then
      q_1 = q(ii,ix5(4,imc))     ! i+1,j+1
      q_2 = q(ii,ix4(4,imc))     ! i+1,j-1
      q_3 = q(ii,ix5(3,imc))     ! i-1,j+1
      q_4 = q(ii,ix4(3,imc))     ! i-1,j-1
   else
      q_1 = -q(1,ix5(4,imc)) -q(4,ix5(4,imc))    ! i+1,j+1
      q_2 = -q(1,ix4(4,imc)) -q(4,ix4(4,imc))    ! i+1,j-1
      q_3 = -q(1,ix5(3,imc)) -q(4,ix5(3,imc))    ! i-1,j+1
      q_4 = -q(1,ix4(3,imc)) -q(4,ix4(3,imc))    ! i-1,j-1
   endif

   partial = idydz *(q_1 - q_2 - q_3  + q_4 )

!   call q_qtensor (q(:,ix5(4,imc)),q_full (:,:,1))  ! j+1,k+1
!   call q_qtensor (q(:,ix4(4,imc)),q_full (:,:,2))  ! j+1,k-1
!   call q_qtensor (q(:,ix5(3,imc)),q_full (:,:,3))  ! j-1,k+1
!   call q_qtensor (q(:,ix4(3,imc)),q_full (:,:,4))  ! j-1,k-1
!   partial = idydz *(q_full (i,j,1) - q_full (i,j,2) - q_full (i,j,3) + q_full (i,j,4))
!idydz*(Qold[i][j+1][k+1][2] - Qold[i][j+1][k-1][2] - Qold[i][j-1][k+1][2] + Qold[i][j-1][k-1][2]) + \
endif




end function partial


!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine surf_moves(q_new)

  use rbf_type
  use variable , only : nn,pp_bulk,Sbulk,chi,gamma_w1,gamma_s,ao_w,ao_s,ix1_surf,normal,nodes,q,lhat,dddx,dddy,dddz,ix2_surf
  use variable , only : dmc
  
  implicit none

  real(dp)    , dimension(3)  :: n
  integer(i4b) :: imc,i,j,k,l,ii
  real(dp)     :: A5,gamma,grad,deltad,surf_L2,partialsurf,anch_deg,delta
  real(dp)    , dimension(3,3)  :: qo_tmp,q_imc,dQx,dQy,dQz,pp,qbar,qdeg,q_tmp_new
  real(dp)    , dimension(3,3)  :: q1x_tmp,q1y_tmp,q1z_tmp
  real(dp)    , dimension(3,3)  :: q2x_tmp,q2y_tmp,q2z_tmp
  real(dp)    , dimension(3,3)  :: q1yn,q1yp,q1xn,q1xp

  logical      :: top
  real(dp)    , dimension(5,nodes)  :: q_new

   A5 = Sbulk/3.d0
ppbulk:if (pp_bulk.lt.1.0)then

   q1x_tmp=0.d0 ; q1y_tmp=0.d0 ; q1z_tmp=0.d0 
   q2x_tmp=0.d0 ; q2y_tmp=0.d0 ; q2z_tmp=0.d0

! Surface Evolution
   do imc =1,nn

 
      if (chi(imc).eq.0.d0)then
         call inv_transformation_tensor( ao_w(:,imc), qo_tmp (:,:))
      else if (chi(imc).eq.1.d0)then
         call inv_transformation_tensor( ao_s(:,imc), qo_tmp (:,:))
      endif

      gamma = gamma_w1*( 1.d0 - chi(imc) )  +  gamma_s*chi(imc) ! value of gamma

      call q_qtensor (q(:,imc),q_imc(:,:))
    
      

      call q_qtensor (q(:,ix1_surf(1,imc)),q1xn (:,:))
      call q_qtensor (q(:,ix1_surf(2,imc)),q1xp (:,:))
      call q_qtensor (q(:,ix1_surf(3,imc)),q1yn (:,:))
      call q_qtensor (q(:,ix1_surf(4,imc)),q1yp (:,:))     



   
      if (ix1_surf(6,imc).ne.0) then         ! Z  + 1

         call q_qtensor (q(:,ix1_surf(6,imc)),q1z_tmp (:,:))
         call q_qtensor (q(:,ix2_surf(6,imc)),q2z_tmp (:,:))

       top =  .false.
    else if (ix1_surf(5,imc).ne.0) then    ! Z - 1

       call q_qtensor (q(:,ix1_surf(5,imc)),q1z_tmp (:,:))
       call q_qtensor (q(:,ix2_surf(5,imc)),q2z_tmp (:,:))
       top = .true.
    endif


! Proyection

    n=normal(3*imc - 2:3*imc)

    
     ! Proyection
    pp(1,1) = 1.d0 - n(1)*n(1); pp(1,2) =      - n(1)*n(2); pp(1,3) =      - n(1)*n(3)
    pp(2,1) =      - n(2)*n(1); pp(2,2) = 1.d0 - n(2)*n(2); pp(2,3) =      - n(2)*n(3)
    pp(3,1) =      - n(3)*n(1); pp(3,2) =      - n(3)*n(2); pp(3,3) = 1.d0 - n(3)*n(3)


    ! Qbar = Q + 1/3*S
    qbar(1,1) = q_imc(1,1) + A5; qbar(1,2) = q_imc(1,2)      ;  qbar(1,3) = q_imc(1,3) 
    qbar(2,1) = q_imc(2,1);      qbar(2,2) = q_imc(2,2) + A5 ;  qbar(2,3) = q_imc(2,3)
    qbar(3,1) = q_imc(3,1);      qbar(3,2) = q_imc(3,2)      ;  qbar(3,3) = q_imc(3,3) + A5

    qdeg = 0.d0


    do i=1,3
       do j=i,3
          
          do k=1,3
             do l=1,3
                qdeg (i,j)= qdeg (i,j)+ pp(i,k)*qbar(k,l)*pp(l,j)       
             enddo
          enddo
       enddo
    enddo

!   gradients

     do i=1,3
        do j=1,3

           
           dQx(i,j) =  (q1xp (i,j) - q1xn (i,j))*dddx
           dQy(i,j) =  (q1yp (i,j) - q1yn (i,j))*dddy

           if(top)then
              dQz(i,j)=dddz * (  3.d0*q_imc(i,j) - 4.d0*q1z_tmp(i,j) + q2z_tmp(i,j) ) 
           else
              dQz(i,j)=dddz * ( -3.d0*q_imc(i,j) + 4.d0*q1z_tmp(i,j) - q2z_tmp(i,j) )              
           endif
        enddo
     enddo


    ii = 0
     do i=1,3
        do j=i,3
           ii = ii +1

           if(i.eq.j) then
              delta = 1.d0/3.d0
           else
              delta =0.d0
           endif

!           dQx(i,j) =
!           dQy(i,j) = 

!           if(top)then
!              dQx(i,j)=dddx * (  3.d0*q_imc(i,j) - 4.d0*q1x_tmp(i,j) + q2x_tmp(i,j) )
!              dQy(i,j)=dddy * (  3.d0*q_imc(i,j) - 4.d0*q1y_tmp(i,j) + q2y_tmp(i,j) )
!              dQz(i,j)=dddz * (  3.d0*q_imc(i,j) - 4.d0*q1z_tmp(i,j) + q2z_tmp(i,j) )
!               
!           else
!              dQx(i,j)=dddx * ( -3.d0*q_imc(i,j) + 4.d0*q1x_tmp(i,j) - q2x_tmp(i,j) )
!              dQy(i,j)=dddy * ( -3.d0*q_imc(i,j) + 4.d0*q1y_tmp(i,j) - q2y_tmp(i,j) )
!              dQz(i,j)=dddz * ( -3.d0*q_imc(i,j) + 4.d0*q1z_tmp(i,j) - q2z_tmp(i,j) )
!               
!           endif


!  L_2 and L_4 terms

           if(ii.eq.1)  surf_L2 = 2.d0 * n(1) * (  dQx(1,1) + dQy(1,2) + dQz(1,3) )
           if(ii.eq.4)  surf_L2 = 2.d0 * n(2) * (  dQx(2,1) + dQy(2,2) + dQz(2,3) )
           if(ii.eq.6)  surf_L2 = 2.d0 * n(1) * (  dQx(3,1) + dQy(3,2) + dQz(3,3) )

           if(ii.eq.2)  surf_L2 = n(1)* (dQx(2,1) + dQy(2,2) + dQz(2,3)  )  + n(2)* (dQx(1,1) +dQy(1,2) + dQz(1,3)  )
           if(ii.eq.3)  surf_L2 = n(1)* (dQx(3,1) + dQy(3,2) + dQz(3,3)  )  + n(3)* (dQx(1,1) +dQy(1,2) + dQz(1,3)  )
           if(ii.eq.5)  surf_L2 = n(2)* (dQx(3,1) + dQy(3,2) + dQz(3,3)  )  + n(3)* (dQx(2,1) +dQy(2,2) + dQz(2,3)  )
          
           partialsurf =  n(1) * (dQx(1,1) + dQy(1,2) + dQz(1,3)) + &
                          n(2) * (dQx(2,1) + dQy(2,2) + dQz(2,3)) + &
                          n(3) * (dQx(3,1) + dQy(3,2) + dQz(3,3))

           grad = dQx(i,j)*n(1) + dQy(i,j)*n(2) + dQz(i,j)*n(3) 


!           if(n==0) surf_L2[0] = 2*nu[i][j][k][0]*(dQx[0]+dQy[1]+dQz[2]);
!           if(n==3) surf_L2[3] = 2*nu[i][j][k][1]*(dQx[1]+dQy[3]+dQz[4]);
!           if(n==5) surf_L2[5] = 2*nu[i][j][k][2]*(dQx[2]+dQy[4]+dQz[5]);
           
!           if(n==1) surf_L2[1] = nu[i][j][k][1]*(dQx[0]+dQy[1]+dQz[2])+nu[i][j][k][0]*(dQx[1]+dQy[3]+dQz[4]);
!           if(n==2) surf_L2[2] = nu[i][j][k][2]*(dQx[0]+dQy[1]+dQz[2])+nu[i][j][k][0]*(dQx[2]+dQy[4]+dQz[5]);
!           if(n==4) surf_L2[4] = nu[i][j][k][2]*(dQx[1]+dQy[3]+dQz[4])+nu[i][j][k][1]*(dQx[2]+dQy[4]+dQz[5]); 
 
!           full = nu[i][j][k][0]*(dQx[0]+dQy[1]+dQz[2]) + \
!                  nu[i][j][k][1]*(dQx[1]+dQy[3]+dQz[4]) + \
!                  nu[i][j][k][2]*(dQx[2]+dQy[4]+dQz[5]);
!           Qnew[i][j][k][n] = Qold[i][j][k][n] - dt*( -1*grad[n]- L2*(0.5*surf_L2[n]-delta3(n)*full) - L3*(0.5*surf_L3[n]-delta3(n)*full) + \ 
!           W*( Qold[i][j][k][n] - Qo[i][j][k][n] ) );           


           anch_deg =0.d0

           do k=1,3
              do l=1,3
                 anch_deg =  anch_deg  + pp(k,l)*qdeg(k,l)
              enddo
           enddo

           if (chi(imc).eq.0.d0)then
              q_tmp_new (i,j) =  q_imc (i,j)  + dmc(1)*( grad +lhat(2)*(0.5d0*surf_L2-delta*partialsurf) - gamma *( (qbar(i,j)-qdeg(i,j)) -deltad(i,j)*(Sbulk- anch_deg)/3.d0 ))
           else if (chi(imc).eq.1.d0)then
              q_tmp_new (i,j) =  q_imc (i,j)  + dmc(1)*( grad +lhat(2)*(0.5d0*surf_L2-delta*partialsurf) - gamma * ( q_imc (i,j) - qo_tmp (i,j) ))
           endif

        enddo
     enddo

 

!call transformation_q( q_tmp_new(:,:), a_new(:,imc) ) ; 
 call qtensor_q (q_tmp_new,q_new(:,imc))
enddo

else
!a_new(:,1:nn) = a (:,1:nn)
q_new(:,1:nn)  = q(:,1:nn)
endif ppbulk



return
 456 format (I6,5f12.8)

end subroutine surf_moves
function dotq(qo,n)
use rbf_type
implicit none
!real(dp) , dimension (3,3):: q
real(dp)  :: dots,dotq
real (dp),dimension(3,3) :: qo
integer :: i,n


dots = 0.d0



if (n.eq.1)then
   do i=1,3
      dots =  dots + qo(1,i)* qo(i,1)
   
   enddo
endif

if (n.eq.2)then
   do i=1,3
      dots =  dots + qo(1,i)* qo(i,2)
   enddo
endif


if (n.eq.3)then
   do i=1,3
      dots =  dots + qo(1,i)* qo(i,3)
   enddo
endif

if (n.eq.4)then
   do i=1,3
      dots =  dots + qo(2,i)* qo(i,2)
   enddo
endif


if (n.eq.5)then
   do i=1,3
      dots =  dots + qo(2,i)* qo(i,3)
   enddo
endif

if (n.eq.6)then
   do i=1,3
      dots =  dots + qo(3,i)* qo(i,3)
   enddo
endif


dotq=dots

end function dotq


function deltad(i,j)
use rbf_type
implicit none
integer i,j
real(dp) :: deltad

if(i.eq.j)then
deltad = 1.d0
else if(i.ne.j)then
deltad=0.d0
endif



end function deltad

!-------------------------------------------------
!-------------------------------------------------!
!-------------------------------------------------
function doubledot(q)
use rbf_type
implicit none
real(dp), dimension(3,3), intent(in) :: q
real(dp) :: doubledot
integer :: i,j 

doubledot = 0.d0
do i=1,3
   do j=1,3
      doubledot = doubledot + q(i,j)*q(j,i)
   enddo
enddo

end function doubledot
